Attachment and detachment rate distributions in deep-bed filtration 
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We study the transport and deposition dynamics of colloids in saturated porous media under un- 
favorable filtering conditions. As an alternative to traditional convection-diffusion or more detailed 
numerical models, we consider a mean-field description in which the attachment and detachment 
processes are characterized by an entire spectrum of rate constants, ranging from shallow traps 
which mostly account for hydrodynamic dispersivity, all the way to the permanent traps associated 
with physical straining. The model has an analytical solution which allows analysis of its properties 
including the long time asymptotic behavior and the profile of the deposition curves. Furthermore, 
the model gives rise to a filtering front whose structure, stability and propagation velocity are ex- 
amined. Based on these results, we propose an experimental protocol to determine the parameters 
of the model. 

PACS numbers: 47.56.+r, 47.55Lm, 47.15.-x 



I. INTRODUCTION 

Many processes in biological systems as well as in the 
chemical and petroleum industry involve the transport 
and nitration of particles in porous media with which 
they interact through various forces jl], ^, |j, These 
interactions often result in particle adsorption and/or 
entrapment by the medium. Examples include filtra- 
tion in the respiratory system, groundwater transport, 
in situ bioremediation, passage of white blood cells in 
brain blood vessels in the presence of jam-1 proteins, 
passage of viral particles in granular media, separation 
of species in chromatography, and gel permeation. The 
particle-medium interactions in these systems are not al- 
ways optimal for particle retention. For example, the 
passage of groundwater through soil often happens un- 
der chemically unfavorable conditions, and as a result 
many captured particles (e.g., viruses and bacteria) may 
be released back to the solution. While filtration un- 
der favorable conditions has been studied and modeled 
extensively^ §, % |, % [ToJ [H], @ g g, we are just 
beginning to understand the process occurring under un- 
favorable conditions. 

Several models have been developed to describe the 
kinetics of particle filtration under unfavorable con- 
ditions. The most commonly used ones are, in 
essence, phenomenological mean-field models based on 
the convection-diffusion equation (CDE) [see Eq. (Q) and 
Sec. 0. Typically, one models the dynamics of free par- 
ticles in terms of the average drift velocity v and the 
hydrodynamic dispersivity A, while the net particle de- 
position rate accounting for particle attachment and 
detachment at trapping sites is a few-parameter function 
of local densities of free and trapped particles. For given 
filtering conditions, the parameters A and v can be de- 
termined from a separate experiment with a tracer, while 
the coefficients of the function can be obtained by fit- 
ting Eq. ([!]) to the breakthrough curves. 



Despite their attractive simplicity, it is widely ac- 
cepted now that the phenomenological models at the 
mean- field level have significant problems. First, the 
depth dependent deposition curves for viruses and bac- 
teria are often much steeper than it would be expected 
if the deposition rates were uniform throughout the 

substrate || [l6[ [T| [01 H H H- This was com- 
monly compensated by introducing the depth-dependent 
deposition rates. The problem was brought to light in 
Ref. p3[ , where it was demonstrated that the steeper- 
than-expected deposition rates under unfavorable filter- 
ing conditions also exist for inert colloids. 

Second, Bradford et al. |24|, pointed out that the 
usual mean-field models based on the CDE, accounting 
for dynamic dispersivity and attachment and detachment 
phenomena, cannot explain the shape of both the break- 
through curves and the subsequent filter flushing. In 
these experiments some particles were retained in the 
medium, and the authors argued for the need to in- 
clude the straining (permanent capture of colloids) in the 
model. Even so, these models may still be insufficient to 
fit the experiments ]2G| ]. 

More elaborate models to describe deep-bed filtration 
have been proposed in Refs. [^7], ||, 29, 3Cj. These mod- 
els go beyond the mean-field description by simulating 
subsequent filter layers as a collection of multiply con- 
nected pipes with a wide distribution of radii, which 
results in a variation in flow speed and also of the at- 
tachment and detachment rates (even straining in some 
cases). The disadvantage of these models is that they 
are essentially computer based: it is difficult to gain an 
understanding of the qualitative properties of the solu- 
tions, without extensive simulations. Furthermore, the 
simulation results suffer from statistical uncertainties. 

In the present work, we develop a minimalist mean- 
field model to investigate filtering under unfavorable con- 
ditions. The model accounts for both a convective flow 
and the primary attachment and detachment processes. 
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Unlike the previous mean-field models of filtration, our 
model contains attachment sites (traps) with different 
detachment rates Bi [see Eq. ([52])], which allows an ac- 
curate modeling of the filtration dynamics over long-time 
periods for a broad range of inlet concentrations. Yet, the 
model admits exact analytical solutions for the profiles of 
the deposition and breakthrough curves which permit us 
to understand qualitatively the effect of the correspond- 
ing parameters and design a protocol for extracting them 
from experiment. 

One of the advantages of our model is that the "shal- 
low" short-lived traps represent the same effect as hydro- 
dynamical dispersivity without generating unphysically 
fast moving particles or requiring an additional boundary 
condition at the inlet of the filter. The "deep" long-lived 
traps allow to correctly simulate long-time asymptotics 
of the released colloids in the effluent during a washout 
stage. The traps with intermediate detachment rates 
determine the most prominent features of breakthrough 
curves. The effect of every trap kind is to decrease the ap- 
parent drift velocity. As attachment and detachment rate 
constants depend on colloid size, we can also account for 
the apparent acceleration of larger particles without any 
microscopic description as in Ref. pl[ . The particle-size 
distribution can be also used to analyze the steeper de- 
position profiles near the inlet of the filter |l6| , 

The paper is organized as follows. In Sec. 
a brief overview of colloid-transport experiments, CDE 
models, and their analytical solutions in simple cases. 
The linearized multiratc convection-only filtration model 
is introduced in Sec. III. The model is characterized by a 
discrete or continuous trap-release-rate distribution; it is 
generally solved in quadratures, and completely in several 
special cases. The results support our argument that 
the hydrodynamic dispersivity can be traded for shallow 
traps. This serves as a basis for the exact solution of 
the full mean-field model for filtration under unfavorable 




conditions introduced in Sec. IV, where we show that a 



large class of such models can be mapped exactly back to 
the linearized ones and analyze their solutions, as well as 
the propagation velocity, structure, and stability of the 
filtering front. We suggest an experimental protocol to 
fit the parameters of the model in Sec. ^ and give our 
conclusions in Sec. VI. 



II. BACKGROUND 

A. Overview of colloid transport experiments 

A typical setup of a colloid-transport experiment is 
shown in Fig. [I]. A cylindrical column packed with sand 
or other filtering material is saturated with water run- 
ning from top to bottom until the single-phase state (no 
trapped air bubbles) is achieved. At the end of this stage, 
colloidal particles are added to the incoming stream of 
water with both the concentration of the suspended par- 
ticles and the flow rate kept constant over time T. This 



is sometimes followed by a filter washout stage in which 
clean water is pumped through the filter. The filtration 
processes are characterized by two relevant experimen- 
tal quantities: the particle breakthrough and deposition 
profile curves. While breakthrough curve represents the 
concentration of effluent particles at the outlet of the 
column as a function of time, deposition curves illustrate 
the depth distribution of concentration of the particles 
retained throughout the column. 



FIG. 1: Schematic of experimental setup in the colloid- 
transport studies. 



B. Convection-diffusion transport model 

As the suspended particles move through the filter- 
ing column, each individual colloid follows its own tra- 
jectory. Consequently, even for small particles that are 
never trapped in the filter, the passage time through the 
column fluctuates. In the case of laminar flows with small 
Reynolds numbers and sufficiently small particles, which 
presumably follow the local velocity lines, the passage 
time scales inversely with the average flow velocity along 
the column v. The effects of the variation between the 
trajectories of particles as well as their speeds can be 
approximated by the velocity-dependent diffusion coeffi- 
cient D = Xv, where A is the hydrodynamic dispersivity 
of the filtering medium. In comparison, the actual dif- 
fusion rate of colloids in experiments is negligibly small. 
Dispersivity is often obtained through tracer experiments 
in which the motion of the particles, i.e., salt ions, which 
move passively through the filter medium without being 
trapped, is traced function of time. 

Overall, the dynamics of the suspended particles along 
the filter can be approximated by the mean- field CDE, 



dC_ 
~dt 



dC 
dx 



Xv 



d 2 C 
dx 2 



-r<i, 



(1) 



where C = C(x,t) is the number of suspended particles 
per unit water volume averaged over the filter cross sec- 
tion at a given distance x from the inlet and Td is the 
deposition rate which may include both attachment and 
detachment processes. 
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C. Issues with the CDE approximation 

The diffusion approximation employed in Eq. (|l|) has 
two drawbacks which could seriously affect the resulting 
calculations if enough care is not used. 

First, while the diffusion approximation works well to 
describe the concentration C(x,t) of suspended parti- 
cles in places where C(x,t) is large, it seems to signif- 
icantly overestimate the number of particles far down- 
stream where C{x, t) is expected to be small or zero. This 
is mainly due to the fact that the diffusion process allows 
for infinitely fast transport, albeit for a vanishingly small 
fraction of particles. In the simple case of tracer dynam- 
ics [rd = in Eq. p )], the general solutions as presented 
in Eqs. (||) and (g) are non-zero even at very large dis- 
tances x — vt ^> 2{Xvt) 1 / 2 . While in many instances this 
may not be crucial, the application of the model to, e.g., 
public health and water safety issues might trigger a false 
alert. 

Second, for the filtering problem one expects the con- 
centration C(x,t) to be continuous, with the concentra- 
tion downstream uniquely determined by that of the up- 
stream. On the other hand, Eq. (|l]) contains second spa- 
tial derivative, which requires in addition to the knowl- 
edge of C(x, t) at the inlet, x = 0, another type of bound- 
ary condition to describe the concentration of particles 
along the column. This additional boundary condition 
could be, e.g., the spatial derivative C'(x,t) at the inlet, 
x = 0, or the outlet, x = L |2|, or the fixed value 
of the concentration at the outlet. We show below that 
fixing a derivative introduces an incontrollable error. On 
the other hand, we cannot introduce a boundary condi- 
tion for the function C(x,t) at the outlet, x = L, as this 
is precisely the quantity of interest to calculate. 

The situation has an analogy in neutron physics 
While neutrons propagate diffusively within a medium, 
they move ballistically in vacuum. A correct calcula- 
tion of the neutron flux requires a detailed simulation of 
the momentum distribution function within a few mean- 
free paths from the surface separating vacuum and the 
medium. In contrast to the filtration theory, for the case 
of neutron scattering, where the neutron distribution is 
stationary it is common to use an approximate boundary 
condition in terms of a "linear extrapolation distance" 
(the inverse logarithmic derivative of neutron density). 

The CDE [see Eq. ([!])] can be solved on a semi-infinite 
interval (x max 3> L) with setting C'(x,t) = at x max 
and calculating the value of C{x,t) at x = L as an ap- 
proximation for the concentration of effluent particles. 
To illustrate this situation, we solve Eq. ([!]) for the case 
of tracer particles, where the deposition rate is set to 
zero, rd = 0. We consider a semi-infinite geometry with 
the initial condition C(x, 0) = and a given concen- 
tration C(0,i) at th e inle t. The corresponding solution 
is presented in Sec. |ll D|. The spatial derivative at the 
boundary given in Eq. (H) is non-zero, time-dependent, 
and rather large at early stages of evolution when the 
diffusive current near the boundary is large. Therefore, 



setting an additional boundary condition for the deriva- 
tive, e.g., C'(0,t) = 0, is unphysical. 

On the other hand, the problem with the boundary 
condition far downstream, C(x max ,i) = 0, x max L, 
can be ill-defined numerically, as this condition is auto- 
matically satisfied to a good accuracy as long as the bulk 
of the colloids has not reached the end of the interval. 



D. Tracer model 

The simplest version of the convection-diffusion equa- 
tion [Eq. ([I])] applies to tracer particles where the depo- 
sition rate is set to zero, = 0, 



dC_ dC _ Xv <PC_ _ 
dt dx dx 2 



(2) 



With the initial conditions, C(x, 0) = 0, the Laplace- 
transformed function C = C(x,p) obeys the equation 



pC + vC' - XvC" = 0, 



(3) 



where primes denote the spatial derivatives, C = 
d x C(x,p). The solution to the above equation is C oc e KX , 
with 



1/1 p 
K±= 2X ± {^ + ^ 



1/2 



(4) 



At semi-infinite interval x > 0, only the solution with 
negative k = K— does not diverge at infinity. Given the 
Laplace-transformed concentration at the inlet, C(0,p), 
we obtain 



C(x,p) = C(0,p) exp I — - x 



J_ p_ 

4A 2 Aw 



1/2N 



(5) 



The inverse Laplace transformation of the above equation 
is a convolution, 



C{x,t) 



dt' C{0,t-t')g(x,t'), 



with the tracer Green's function (GF) 



g{x,t) 



1 



2( 7 rAv) 1 /2 £3/2 



exp 



{x - vt)' 
AXvt 



(G) 



(7) 



In the special case C(0,<) — Cq =const, the integration 
results 



C = 9l U + erf f tV ~ X I + erfc [ tv + x 

(8) 

where erfc(z) = 1 — erf(z) is the complementary error 
function. 

We note that the spatial derivative of the solution of 
Eq. m) at x = is different from zero. Indeed, it depends 
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on time and is divergent at small t, implying an unphys- 
ically large diffusive component of the particle current, 



C(0,t) = ^ 



erfc (a) 



2A (nivX) 1 / 2 



a 



4A- (9) 



In the presence of the straining term, = A N C in 
Eq. (g), the GF can be obtained from Eq. (Q) by intro- 
ducing exponential decay with the rate AqNq, 



SO: t) 



-A N a t 



exp 



{x - vtf 
AXvt 



(10) 



Note that we wrote the straining rate as a product of the 
capture rate Aq by infinite-capacity "permanent" traps 
with the concentration Nq per unit volume of water. Such 
a factorization is convenient for the non-linear model pre- 
sented later in Sec. IV. The same notations are employed 



throughout this work for consistency. 



III. LINEARIZED MEAN-FIELD FILTRATION 
MODEL 

In this section we discuss the linearized convection- 
only multitrap filtration model, a variant of the multi- 
rate CDE model first proposed in Ref. |54|. Our model 
is characterized by a (possibly continuous) density of 
traps as a function of detachment rate [see Eq. (p3|)]. 
Generically, continuous trap distribution leads to non- 
exponential (e.g., power-law) asymptotic forms of the 
concentration in the effluent on the washout stage. 

The main purpose of this section is to demonstrate 
that "shallow" traps with large detachment rates have 
the same effect as the hydrodynamic dispcrsivity in CDE. 
In addition, the obtained exact solutions will be used 



in Sec. |IVj as a basis for the analysis of the full non- 
linear mean-field model for filtration under unfavorable 
conditions. 



with n\ = ni(x, t) as the auxiliary variable describing the 
average number of particles in a trap, N\ as the number 
of traps per unit water volume, A\ as the trapping rate, 
and B\ as the release rate. The particular normalization 
of the coefficients is chosen to simplify the formulation of 
models with traps subject to saturation in Sec. IV. 

To simulate dispcrsivity where all time scales are in- 
versely proportional to propagation velocity, we must 
choose both A\ and B\ proportional to v. The corre- 
sponding parameter a in A\ = av has a dimension of 
area and can be viewed as a trapping cross section. The 
length I in the release rate B\ = v/£ can be viewed as 
a characteristic size of a stagnation region. On general 
grounds we expect a oc £ 2 with i on the order of the grain 
size. 



B. Single-trap model with straining. 

To illustrate how shallow traps can provide for dis- 
persivity in convection-only models, let us construct the 
exact solution of Eq. (pd|). In fact, it is convenient to 
consider a slightly generalized model with the addition 
of straining, 



C + vC + N±ni 



-A N C, fn = A 1 C-B 1 m. (12) 



With zero initial conditions the Laplace transformation 
gives for C = C(x,p), 



p + A N + ) C + vC' = 0. (13) 



The boundary value for Laplace-transformed C(x, t) at 
the inlet is given by C(Q,p). With initially clean filter, 
C(x,0) = n(x,0) = 0, and a given free particle con- 
centration C(0, t) at the inlet, the solution to the linear 
one-trap convection-only model with straining [Eq. (p^)] 
is a convolution of the form presented in Eq. (|g) with the 
following GF §5): 



A. Shallow traps as a substitute for diffusion 

To rectify the problems with the diffusion approxi- 
mation noted previously, we suggest an alternative ap- 
proach for the propagation of particles through the fil- 
tering medium. Instead of considering the drift with an 
average velocity with symmetric diffusion-like deviations 
accounting for dispersion of individual trajectories, we 
consider the convective motion with the maximum ve- 
locity v. The random twists and turns delaying the indi- 
vidual trajectories are accounted for by introducing Pois- 
sonian traps which slow down the passage of the majority 
of the particles through the column. In the simplest case 
suitable for tracer particles, the relevant kinetic equations 
read as follows: 



C + vC + N 1 h 1 = Q, hi = A 1 C-B 1 m, 



(11) 



g(x,t) = e-^ v - B ^- x ^h(t-x/v) 

+6(t - x/v)— Ji(Ct) 



(14) 



where f3 = A N + A1N1 is the clean-bed trapping rate, 
9{z) is the Heaviside step- function, and h((t) is the mod- 
ified Bessel function of the first kind with the argument 



[A^xB^tv ~ x)x] 1/2 . 



(15) 



The singular term with the 8 function S(t — x/v) in 
Eq. (|l4|) represents the particles at the leading edge 
which propagate freely with the maximum velocity v 
without ever being trapped. The corresponding weight 
cxp(— (3x/v) decreases exponentially with the distance 
from the origin. 
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Sufficiently far from both the origin and from the lead- 
ing edge, where the argument (t [Eq. (|l5|)1 of the Bessel 
function is large, we can use the asymptotic form, 



hit) 



Subsequently, Eq. (h4h becomes 



(2ttC) 1 /2 



; c [1 + OiC 1 )] ,RcC>0. (16) 



g(x, t) 



-A N x/v 



Bit 



1/4 



27r l/2 r 3/4 



cxp-(^-^) 2 , (17) 



where r = Bi(t — x/v) is the dimensionless retarded time 
in units of the release rate, and £ = A\N\x/v is the di- 
mensionless distance from the origin in units of the trap- 
ping mean free path. 

The correspondence with the GF in Eq. ( [To| ) for the 
CDE with linear straining [or Eq. (0) for the CDE tracer 
model in the case of no permanent traps, Nq = 0] can be 
recovered from Eq. ( |l7| ) by expanding the square roots 
in the exponent around its maximum at £ = r, or x — 
vot, with the effective velocity vq = vB\/{B\ + NiAi). 
Specifically, suppressing the prefactor due to straining, 
[No = in Eq. (p"7|)], we obtain for the asymptotic form 
of the exponent at large t, 



l(x, t) oc exp • 



(x - V t) 2 



with the effective dispersivity coefficient [cf. Eq. 



An 



(N 1 A 1 +B 1 ) 



(18) 



(19) 



The approximation is expected to be good as long as 
both x and t are large compared to the width of the bell- 
shaped maximum. 

The actual shapes of the corresponding GFs, Eqs. (0) 
and ( pT[ ) in the absence of permanent traps, iVn = 0, 
are compared in Fig. ||. While the shape differences are 
substantial at small t, they disappear almost entirely at 
later times. 




FIG. 2: (Color online) Comparison of the spatial depen- 
dence of the GFs for the tracer model implemented as the 
convection-diffusion equation [Eq. (|l|)] with r& = (solid 
lines) and the single-trap convection model [Eq. ([n])] (dashed 
lines). Specifically, we plot Eq. (Q) and the regular part of 
Eq. ( |l4j ) with iVo = 0, using identical values of v = Vo = 1 
and A = Ao = 1 and the release rate Bi = 1/2 (half the max- 
imum value at these parameters) at t = 2, 4, 8, 16, 32. Once 
the maximum is sufficiently far from the origin, the two GFs 
are virtually identical (see Sec. IIIB). 



The corresponding solution can be obtained in quadra- 
tures in terms of the Laplace transformation. With the 
initial condition, C(x, 0) = rii(x, 0) = and a given time- 
dependent concentration at the inlet, C(0, t) — Cn(i), the 
result for C(x,t) is a convolution of the form presented in 
Eq. (^) with the GF given by the inverse Laplace trans- 
formation formula, 



g{x,t) 



dp p [t — x/v — xT,(p)/v] 



2iri 



with the response function 



dB p(B) 
P + B 



(21) 



(22) 



C. Multitrap convection-only model 

Even though the solutions of the single-trap model cor- 
respond to those of the CDE [Eq. (^)], the model pre- 
sented in Eq. (|lj) is clearly too simple to accurately de- 
scribe filtration under conditions where trapped particles 
can be subsequently released. At the very least, in addi- 
tion to straining and "shallow" traps that account for the 
dispersivity, describing the experiments J24|, |25| requires 
another set of "deeper" traps with a smaller release rate. 

More generally, consider a linear model with m types 
of traps differing by the rate coefficients Ai, Bi, 

m 

C + vC + ^Nihi = 0, n, = AiC-Bim. (20) 

t=i 



Here we introduced the effective density of traps, 

m 

p{B)=Y J A t N l 5{B-B i ), (23) 

i=l 

corresponding to various release rates. 

The general structure of the concentration profile can 
be read off directly from Eq. (pl|). It gives zero for 
t < x/v, consistent with the fact that v is the maximum 
propagation velocity in Eq. (p0|). The structure of the 
leading-edge singularity (the amplitude of the 5 function 
due to particles which never got trapped) is determined 
by the large-p asymptotics of the integrand in Eq. ([2l]). 
Specifically, GF ([H]) can be written as 

g(x, t) = e- f3x / v 6(t - x/v) + 9(t - x/v)g ieg (x, t), (24) 
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where (3 = lim^oc p£(p) = £\ [cf. Eq. (jlj)] is the 

clean-bed trapping rate, and (? reg is the non-singular part 
of the GF. 

Similarly, the structure of the diffusion-like peak of 
the GF away from both the origin and the leading edge 
is determined by the saddle point of the integrand in 
Eq. (|2l]) at small p. Assuming the expansion S(p) = 
E(0) — Sip + 0(p 2 ) and evaluating the resulting Gaus- 
sian integral around the saddle point at 



we obtain 



g(x,t) 



t — x/Vq 

2xT,i/v ' 



2(ttJ: 1 x/v) 1 / 2 



1 + S(0)' 



-(t-x/v )'/(4Eix/v) 



(25) 



(26) 



The exponent near the maximum can be approximately 
rewritten in the form of that in Eq. (Q) , with the effective 
dispcrsivity 



Ao — — 2ji 

V 



[1 + E(0)] S 



(27) 



For the case of one trap, m = 1, the expressions for 
the effective parameters clearly correspond to our ear- 
lier results of Eqs. (|§) and @. Note that the precise 
structure of the exponent and the prefactor in Eq. (|2?j ) 
is different from those in Eq. (|l^) which was obtained by 
a more accurate calculation. 

The effective diffusion approximation [Eq. (^6|)] is ac- 
curate for large x near the maximum as long as the inte- 
gral in Eq. (Ell) remains dominated by the saddle-point 
in Eq. (25). In particular, the poles of response function 
( p2| ) must be far from p*. This is easily satisfied in the 
case of "shallow" traps with large release rates Bi 3> \p*\. 

On the other hand, this condition could be simply vio- 
lated in the presence of "deep" traps with relatively small 
Bi. Over small time intervals compared to the typical 
dwell time B~ , these traps may work in the straining 
regime in which they would not contribute to the effec- 
tive dispersivity. This situation may be manifested as an 
apparent time-dependence of the effective drift velocity 
Vq and/or the dispersivity Ao- 



D. Model with a continuous trap distribution 



The multitrap generalization given in Eq. (|20|) for fil- 
tration is clearly a step in the right direction if we want 
an accurate description of the filtering experiments. 

Indeed, apart from the special case of a regular array 
of identical densely-packed spheres with highly polished 
surfaces, one expects the trapping sites (e.g., the contact 
points of neighboring grains) to differ. For small parti- 
cles such as viruses, even a relatively small variation in 
trapping energy could result in a wide range of release 
rates Bi differing by many orders of magnitude pq, E9|. 



Under such circumstances, it is appropriate to consider 
mean-field models with continuous trap distributions. 

Here we only consider a special case of a continuous 
distribution of the trap parameters, Ai and Bi, such that 
the release-rate density in Eq. (E2) has an inverse-square- 
root singularity, p{B) = Pi^/T^B 1 / 2 ), with the release 
rates ranging from infinity all the way to zero. The cor- 
responding response function (p2|) could be expressed as 



S(P) = Pl/2/p 1/2 - 



(28) 



The inverse Laplace transform [Eq. (pl|)] gives the follow- 
ing GF: 



g(x,t) 



XPl/2 



2 p? /2 /(4- 2 



>0(T), 



t . 



2y / 7TUT 3 / 2 V 

_ (29) 

Note that, in accordance with Eq. (J24|), there is no 
leading-edge S function near t = x/v as the expression 
for the corresponding trapping rate (3 diverges. Because 
of the singular behavior of S(p) at p = 0, there is no 
saddle-point expansion of the form given in Eq. (|25|). 
Thus, there is no Gaussian representation analogous to 
Eq. (|2^): at large t, the maximum of the GF is located 
a t x max = vp 1 / 2 i2t) 1 / 2 , which is also of the order of the 
width of the Gaussian maximum. The GF [Eq. Q29)] for 
two representative values of p\ji is plotted in Fig. |3|. 




FIG. 3: (Color online) Spatial dependence of the GF [Eq. (g9|)] 
for the model presented in Eq. Q20J) with continuous distribu- 
tion of trap parameters corresponding to inverse-square-root 
singularity in the response function [see Eqs. ( |22| and (p8[))]. 
Dashed lines show the GF at p 1 i 2 = 0.25, while solid lines 
present the same GF at pi/ 2 — 1 multiplied by the factor of 
8. We chose t = 2,4, 8, 12 as indicated in the plot. Unlike in 
Fig. ^, due to abundance of traps with long release time, the 
GFs do not asymptotically converge toward a Gaussian form. 



We also note that for large t at any given x, Eq. (|2£ 
has a power-law tail oc i~ 3//2 . This property is generic for 
continuous trap density distributions leading to small-p 
power-law singularities in For example, taking the 

density of the release rates as a power law in B, 



P(B) 



sin(7rs) p s 
n ~B~ S 



(30) 
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where s is the corresponding exponent. < s < 1, we 
obtain = p s p s , and the large-t asymptotic of the 
GF at a fixed finite x scales as 



ticlcs with concentration C, 



g(x, t) cx t s 



(31) 



Such a power law is an essential feature of continuous 
distribution (^0|) of the detachments rates; it cannot be 
reproduced by a discrete set of rates Bi which always 
produce an exponential tail. 



IV. FILTRATION UNDER UNFAVORABLE 
CONDITIONS 

A. Multitrap model with saturation 

The considered linearized filtration model presented by 
Eq. (^0|) can be used to analyze filtration of identical 
particles in small concentrations and over limited time 
interval as long as the trapped particles do not affect the 
filter performance. However, unless the model is used 
to simulate tracer particle dynamics in which no actual 
trapping occurs, it is unlikely that the model remains 
valid as the number of trapped particles grows. 

Indeed, one expects that a trapped particle changes 
substantially the probability for subsequent particles to 
be trapped in its vicinity. Under favorable filtering condi- 
tions characterized by filter ripening |36L |37j , the proba- 
bility of subsequent particle trapping increases with time 
as the number of trapped particles ni grows. On the 
other hand, under unfavorable filtering conditions, where 
the Debye screening length is large compared to the trap 
size £, for charged particles one expects trapping proba- 
bilities Ai{ni) to decrease with rn. 

If repulsive force between particles is large, we can as- 
sume that only one particle is allowed to be captured in 
each trap. Subsequently, a single trap can be character- 
ized by an attachment rate Ai when it is empty and a de- 
tachment rate Bi when it is occupied, and the mean-field 
trapping/release dynamics for a given group of trapping 
sites can be written as 



= CAi(l -m) -Bim 



(32) 



Note that this equation is non-linear because it contains 
the product of CVij. 

Previously, similar filtering dynamics was considered 
in a number of publications (see Refs. [^4| and [^5) and 
references therein). In the present work, we allow for 
a possibility of groups of traps differing by the rate pa- 
rameters Ai and Bi. The distribution of rate parameters 
can also be viewed as an analytical alternative of the 
computer-based models describing a network of pores of 
varying diameter |7[ |9],[30). 

Our mean-field transport model is completed by 
adding the kinetic equation for the motion of free par- 



C + vC + ]T Nihi = 0, 



(33) 



which has the same form as the linearized equations 
[Eq. @] considered in Sec. |HID . 

We note that for shallow traps with large release rates 
Bi , the non-linearity inherent in Eq. ([32]) is not important 
for sufficiently small suspended particle concentrations 
C. Indeed, if C is independent of time, the solution of 
Eq. (B3) saturates at 



i{C) = 



B l + CA, 



(34) 



For small free-particle concentration C, or for any C and 
large enough Bi, the trap population is small compared 
to 1, and the non- linear term in Eq. (|3^) can be ignored. 

Therefore, as discussed in relation with the linearized 
multitrap model [sec Sec. Ill A and Eq. (p0|)], the effect 



of shallow traps is to introduce dispersivity of the arrival 
times of the particles on different trajectories. For this 
reason, we are free to drop the dispersivity term [cf. the 
CDE model, Eq. (0)], and use a simpler convection-only 
model ([}j|) with several groups of traps with density Ni 
per unit water volume, characterized by the relaxation 
parameters Ai and Bi. 



B. General properties: Stable filtering front 



The constructed non-linear equations [Eqs. ([32|) and 
(|33|)] describe complicated dynamics which is difficult to 
understand in general. Here, we introduce the front ve- 
locity, a parameter that characterizes the speed of dete- 
rioration of the filtering capacity. 

Consider a semi-infinite filter, with the filtering 
medium initially clean, and the concentration C(0,t) = 
Ca of suspended particles at the inlet constant. Af- 
ter some time, the concentration of deposited particles 
near the inlet reaches the dynamical equilibrium u^Ca) 
[Eq. (HJ)] and, on average, the particles will no longer be 
deposited there. At a given inlet concentration, the fil- 
tering medium near the inlet is saturated with deposited 
particles. On the other hand, sufficiently far from the 
inlet, the filter is still clean. On general grounds, there 
should be some crossover between these two regions. 

The size of the saturated region grows with time [see 
Fig. The corresponding front velocity va = v(Ca) can 
be easily calculated from the particle balance equation, 



v A C A + v A ^ N i n i{CA) = vC A - 



(35) 



This equation balances the number of additional particles 
needed to increase the saturated region by Sx = VA$t on 
the left, with the number of particles brought from the 
inlet on the right [see Fig. |4| . The same equation can also 
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FIG. 4: Solid line shows the free particle concentration near a 
filtering front. Dashed line shows the front shifted by Ax; the 
additional free and trapped particles in the shaded region are 
brought from the inlet [see Eq. (|35|)], See Eq. (j5^) for exact 
front shape. 



be derived if we set C = C(x — v^t), rij = n^x — v^t) and 
integrate Eq. (^3|) over the entire crossover region. The 
trapped particle density saturates as given by Eq. (34), 
and the resulting front velocity is 



FIG. 5: Free particle concentration C(x,t) with two filtering 
fronts. The initial front moves on the background of clean fil- 
ter and leaves behind the equilibrium filtering medium with 
C = Ca- The secondary front with higher inlet concentration 
C'b is moving on partially saturated medium. With nonlin- 
earity as in Eq. (p2|), the secondary front is always faster, 
vab > va; the two fronts will eventually coalesce into a single 
front. 



C. Exactly solvable case 

1. General solution 



v{C A ) 



(36) 



A,C A + Bi 



This is a monotonously increasing function of Ca'- larger 
inlet concentration Ca leads to higher front velocity, 
which implies that the filtering front is stable with re- 
spect to perturbations. Indeed, in Appendix we show 
that the velocity vab of a secondary filtering front with 
the inlet concentration Cb > Ca (see Fig. ||), moving 
on the background of equilibrium concentration of free 
particles Ca, is higher than va, i.e., vab > va- Thus, 
if for some reason the original filtering front is split into 
two parts, moving with the velocities v A and vab, the 
secondary front will eventually catch up, restoring the 
overall front shape. 

We emphasize that the existence of the stable filter- 
ing front is in sharp contrast with the linearized filtering 
problem [see Eq. (p0|)] , where the propagation velocity 
vq [Eq. (25)] is independent of the inlet concentration, 
and any structure is eventually washed out dispersively 
(the width of long-time GF does not saturate with time). 
Also, in the case of the filter ripening, the nonlinear term 
in Eq. ([32j) will be negative and thus would prohibit the 
filtering front solutions due to the fact that the secondary 
fronts move slower, vab < va- The non-linear problem 
with saturation is thus somewhat analogous to Korteweg- 
de Vries solitons where the dispersion and nonlincarity 
compete to stabilize the profile fl38[ |39[ . 



Compared to the linear case presented in Sec. [II, the 
physics behind the non-linear equations [Eqs. (|3 2|) and 
(p3|)] is much more complicated. However, the struc- 
ture of these equations immediately indicates that non- 
lincarity reduces filtering capacity because trapping sites 
could saturate in this model [see Eq. (S3)]. While the 
relevant equations can also be solved numerically, a thor- 
ough understanding of the filtering system, especially 
with large or infinite number of traps, is difficult to 
achieve. 

To gain some insight about the role of the different 
parameters in the filtering process, we specifically focus 
on the non-linear models presented by Eqs. ( |32| ) and ( |3"3| ) 
which can be rendered into a linear set of equations, very 
similar to the linear multitrap model [Eq. (p0|)]. To this 
end, we consider the case where all trapping sites have 
the same trapping cross sections, that is, all A; = A in 
Eq. (|32|). If we introduce the time integral 



u(x, t) 



C(x,t')dt', 



(37) 



then Eq. (|3 
written as 



after a multiplication by expAu can be 



Ot (e Au ) 



d t {n t e Au ) + Bi (me Au 
Clearly, these are a set of linear equations 

in + Bidi = w, 
with the following variables 



w = w(x, t) = e 



Au 



di(x, t) = Hi W. 



(38) 
(39) 
(40) 
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Note that Eq. (|33|) can also be written as a set of lin- 
ear equations in terms of these variables. If we inte- 
grate Eq. over time, we find 



it + vu + NiTii = 0, 



(41) 



where we assumed initially clean filter, C(x,0) = 
rii(x,0) = 0. Considering that w = Aiiw and w' = 
Au' w, we obtain 



w + vw 



' + A J2 Nidi = 0. 



(42) 



The main difference of the linear Eqs. ( |39| ) and ( |42] ) 
from Eqs. ( p(i|) is in their initial and boundary conditions, 



w(x, 0) 
w(Q,t) 



1, a,i(x, 0) = 0, 



„Au (t) 



U (t) 



(43) 

dt'C(0,t'). (44) 



Note that with the time-independent concentration of the 
particles in suspension at the inlet, i.e., C(0,t) = Co, 
boundary condition (H) gives a growing exponent, 



(45) 



w {t) = w(0,t) = e AOot . 



The derived equations can be solved with the use of 
the Laplace transformation. Denoting w = w(x,p) = 
C p {w(t)} and eliminating the Laplace-transformed trap 
populations fii(x,p) = £ p {n,;(x, t)}, we obtain 



(pw - 1) [1 + E(p)] + vw' = 0, Z(p) = Aj2 



N 4 . 



P + Bt 
(46) 

The response function is identical to that in 

Eq. (p2]), and for the case of continuous trap distribu- 
tion we can also introduce the effective density of traps, 
p{B) = AY,i Ni5(B - Bi). The solution of Eq. @ and 
the Laplace-transformed boundary condition [Eq. (44)] 
becomes 



1 



w 



w Q (p) 



1 



-[l+T,(p)]px/v 



(47) 



where w(0,p) = wq(p). Employing the same notation as 
in Eq. (|2l|), the real-time solution of Eqs. (|39| ) and j42] ) 
with the boundary conditions [Eqs. ([13]) and (Q)] can 
be written in quadratures, 

w(x,t) = l+[ dt' [w (t-t')-l] g(x,t'). (48) 



The time-dependent concentration can be restored from 
here with the help of logarithmic derivative, 



C{x,t) 



1 d \nw(x, t) 



A 



Of 



(49) 



2. Structure of the filtering front 

In the special case C(0, t) = Cq = const, the inte- 
grated concentration [Eq. (|37|)1 is linear in time at the 
inlet, uo(t) = Cot, and w(0,t) grows exponentially [see 
Eq. (p5])]. This exponent determines the main contribu- 
tion to the integral in Eq. ( f48| ) for large t and x. In- 
deed, in this case wc can rewrite Eq. ( f48| ) exactly as 
w(x,t) = 1 + J (C ) - J(0), where 



J (C ) 



= e ACot 



-AC of 



g(x,t'). 



(50) 



Note that J(0) is proportional to the solution of the lin- 
earized equations [Eq. (20)] with time- independent inlet 
concentration C(0,t) = const [see Eq. (§)]. The corre- 
sponding front is moving with the velocity vp [Eq. (|2^)] 
and is widening over time [Eqs. (§|) and ©]. Thus, for 
x/uo — t positive and sufficiently large, this contribution 
to w(x, t) is small and can be ignored. In the opposite 
limit of large negative x/vq — t, J(0) = 1, which exactly 
cancels the first term in Eq. (|48|). 

On the other hand, the term J(Cq) grows exponen- 
tially large with time. At large enough t, the integration 
limit can be extended to infinity, and the integration in 
Eq. (|50| ) becomes a Laplace transformation, thus 



w(x, t) 



s AC t 



dt' e- ACot ' g(x,t') 



e Pat e -[l+Y,(p a )}p Q x / v 



p = AC . (51) 



This results in the following free-particle concentration 
[see Eq. ©], 



C(x,t) 



C 



e [x/v(C )-t]ACa _|_ i ' 

and the occupation of the ith trap [Eqs. (^ 

.4 



rii(x,t) = 
with the front velocity 
v(C ) 



B, + ACo 



C(x,t), 



1 + S(AC )' 



(52) 
and ©], 
(53) 

(54) 



Note that this coincides exactly with the general case 
presented in Eq. (|3q) if we set all A4 — A. 



3. Filtering front formation 

The approximation in Eq. ( |5l| ) is valid in the vicinity of 
the front, \x/v(Co) — 1\ < (ACo)^ 1 , as long as x/vo — t is 
positive and large. Since v(Co) > vq = f(0), this implies 



1 

Vo 



1 

AC~o 



» 



^C : 



(55) 
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which provides an estimate of the distance from the 
outlet where the front structure [Eqs. (|||) and (^||)] is 
formed. The exactness of the obtained asymptotic front 
structure can be verified directly by substituting the ob- 
tained profiles in Eqs. ( |3^ ) and ([33). 

The exact expressions in Eqs. (4q) and ( f49| ) for the 
free-particle concentration can be integrated completely 
in some special cases. Here we list two such results and 
demonstrate the presence of striking similarities in the 
profiles C(x,t) between different models, despite their 
very different rate distributions. Furthermore, we show 
that the corresponding exact solutions [Eq. converge 
rapidly toward the general filtering front [Eq. ( |52|)]. 
Single-trap model with straining. In Sec. [MB , we 
found the explicit expression [Eq. ([Ly)] for the GF in 
the case of the linear model for two types of trapping 
sites with rates A\ and B\ and permanent sites with the 
capture rate A . The resulting GF (with Ai = A = A 
and B\ = B) can be used in Eq. ( |48| ) to construct the 
solution for the corresponding model with saturation, 



C- 
n 



vC + N h 
= AC(l - n ) 



Nxfl! 

hi = 



= 0, 
ACQ. 



(56) 
Bim. (57) 



Let us consider the special case of the inlet concentra- 
tion, C(0,<) — Cq6(T — t)6(t), constant over the interval 
< t < T, and zero afterwards. The function wo(t) [see 
Eq. @)] is, then 



Wo(t) = expL4Co min(£, T)], 
and the integration in Eq. (jig) gives 
w = 1 + e~ K [W(t) - e ACoT W(t - T)] , 

w(t) = e(t-o{ 



(58) 



(59) 



,AC (t-f) 



dTe -B l{ r-0 



,AC {t~ 



fo(Cr)}, (60) 



where £ = x/v and £ r is given in Eq. (|l5|). The con- 
centration of free particles, C{x, t), can be now obtained 
through Eq. (^9|). The step function 0(t — x/v) included 
in w indicates that it takes at least t = x/v for a particle 
to travel a distance x. 

Figure ^| illustrates C(x, t) as a function of distance, x, 
at a set of discrete values of time t = 1, 2, ...,16. The 
model parameters as indicated in the caption were ob- 
tained by fitting the response function S(p) = AN /p + 
ANi/ip + Bx) at the interval 0.5 < p < 5.0 to that of the 
model with the continuous trap distribution (sec Fig. |^). 
The solid lines show the curves for t < T, while the 
dashed lines correspond to t > T; they have a drop in the 
concentration near the origin consistent with the bound- 
ary condition at the inlet. The exact profiles show excel- 
lent convergence toward the corresponding front profiles 
computed using Eq. ( |52| ) (symbols). 
Model with square-root singularity. Let us now con- 
sider the non-linear model, [Eqs. (53|) and (p2|)] with the 
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FIG. 6: (Color online) Formation of the filtering front for 
the single-trap filtering model with straining [Eq. (pq)]. Lines 
show the free-particle concentration C(x,t) extracted from 
Eq. @ with T = 10, A = v = Co = 1, No = 0.388, Ni = 
3.60, and B 1 = 4.97, for t = 1, 2, ... ,16. Symbols show the 
front solution [Eq. for t > 10 with the front velocity 

[Eq. (§)]. 



inverse-square-root continuous trap distribution, produc- 
ing the response function given in Eq. (p8|). The model is 
exactly solvable if we set all Aj = A, while allowing the 
trap densities Ni vary with B appropriately. 

The solution for the auxiliary function w correspond- 
ing to the inlet concentration C(0,i) constant on an in- 
terval of duration T is obtained by combining Eqs. (|4g| ) 
and (H), with the relevant GF [Eq. (§9|)]. The result- 
ing ^-dependent curves C(x, t) at a set of discrete time 
values arc shown in Fig. (R), along with the correspond- 
ing asymptotic front profiles (symbols), for a parameter 
set as indicated in the caption. The solid lines show the 
curves for t < T. The dashed lines are for t > T; they 
display a drop of the concentration near the origin con- 
sistent with the boundary condition at the inlet. Again, 
the time-dependent profiles show gradual convergence to- 
ward front solution J52|). 

Note that the profiles in Figs. ^| and [7] are very simi- 
lar even though the corresponding trap distributions dif- 
fer dramatically. This illustrates that parameter fitting 
from a limited set of breakthrough curves is a problem 
ill-defined mathematically. The complexity and ambigu- 
ity of the problem grow with increasing number of traps. 
In Sec. we suggest an alternative computationally sim- 
ple procedure for parameter fitting using the data from 
several breakthrough curves differing by the input con- 
centrations. 



V. EXPERIMENTAL IMPLICATIONS 

The suggested class of mean-field models is character- 
ized by a large number of parameters. In the discrete 
case, these are the trap rate constants Aj, Bi and the cor- 
responding concentrations Ni along with the flow velocity 
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FIG. 7i_f Color online) As in Fig. g but for filtering model (g3|) , 
Eq. ( p2|) with continuous inverse-square-root trap distribu- 
tion [Eq. (psl)]. Parameters are A = v = C n = P1/2 = 1, 
T = 10. Symbols show the front solution [Eq. (J52[)] for t > 10 
with front velocity (|54[). The raising parts of the curves are 
almost identical with those in Fig. H, while there are some 
quantitative differences in the tails, consistent with the ex- 
ponential vs power-law long-time asymptotics of the corre- 
sponding solutions. 



v. In the continuous case, the filtering medium is char- 
acterized by the response function £(p) [see Eq. (p2)]- In 
our experience, two or three sets of traps are usually suffi- 
cient to produce an excellent fit for a typical experimental 
breakthrough curve (not shown). This is not surprising, 
given the number of adjustable parameters. On the other 
hand, from Eq. (54) it is also clear that the obtained pa- 
rameters would likely prove inadequate if we change the 
inlet concentration. The long-time asymptotic form of 
the effluent during the washout stage would also likely 
be off. 



One alternative to a direct non-linear fitting is to use 
our result given in Eq. (Q) [or Eq. (p6|)] for the filtering 
front velocity as a function of the inlet concentration, 
Co- With a relatively mild assumption that all trapping 
rates coincide, Ai = A, one obtains the entire shape of 
the filtering front [Eq. (p2[)]. Thus, fitting the front pro- 
files at different inlet concentrations Co to determine the 
parameter A and the front velocity v(Cq) can be used to 
directly measure the response function S(p). 

The suggested experimental procedure can be summa- 
rized as follows, (i) One should use as long filtering 
columns as practically possible in order to achieve the 
front formation for a wider range of inlet concentrations, 
(ii) A set of breakthrough curves C(L, t) for several con- 
centrations Co at the inlet should be taken, (iii) For 
each curve, the front formation and the applicability of 
the simplified model with all Ai = A should be verified 
by fitting with the front profile [Eq. (p^)]. Given the 
column length, each fit would result in the front veloc- 
ity i>(Co), as well as the inverse front width p = ACq. 
(iv) The resulting data points should be used to recover 
the functional form of T,(p) and the solution for the full 



model. 

It is important to emphasize that the applicability of 
the model can be controlled at essentially every step. 
First, the time-dependence of each curve should fit well 
with Eq. (52). Second, the values of the trapping rate A 
obtained from different curves should be close. Third, the 
computed washout curves should be compared with the 
experimentally obtained breakthrough curves. The ob- 
tained parameters, especially the details of £(p) for small 
p, can be further verified by repeating the experiments 
on a shorter filtering column with the same medium. 



VI. CONCLUSIONS 

In this paper, we presented a mean-field model to in- 
vestigate the transport of colloids in porous media. The 
model corresponds to the filtration under unfavorable 
conditions, where trapped particles tend to reduce the fil- 
tering capacity, and can also be released back to the flow. 
The situation should be contrasted with favorable filter- 
ing conditions characterized by filter ripening. These two 
different regimes can be achieved, e.g., by changing pH 
of the media if the colloids are charged. The unfavorable 
filtering conditions are typical for filtering encountered in 
natural environment, e.g., ground water with biologically 
active colloids such as viruses or bacteria. 

The advantages of the model are twofold. It not only 
fixes some technical problems inherent in the mean-field 
models based on the CDE but also admits analytical so- 
lutions with many groups of traps or even with a contin- 
uous distribution of detachment rates. It is the existence 
of such analytical solutions that allowed us to formulate 
a well-defined procedure for fitting the coefficients. Ulti- 
mately, this improves predictive capability and accuracy 
of the model. 

The need for the attachment and detachment rate dis- 
tributions under unfavorable filtering conditions has al- 
ready been recognized in the field [^J, ^5], ^6|. Previously 
it has been implemented in computer-based models in 
terms of ad hoc distributions of the pore radii [^7| ^{], po| . 
Such models could result in good fits to the experimental 
breakthrough curves. However, we showed in Sec.[y] that 
the relevant experimental curves are often insensitive to 
the details of the trap parameter distributions, especially 
on the early stages of filtering. 

On the other hand, our analysis of the filtering front 
reveals that the front velocity as a function of the inlet 
colloid concentration, v(Cq) [Eq. (|36|)], is primarily de- 
termined by the distribution of the attachment and de- 
tachment rates characterizing the filtering medium. We, 
indeed, suggest that the filtering front velocity is one of 
the most important characteristics of the deep-bed filtra- 
tion as it is directly related to the loss of filtering capacity. 

We have developed a detailed protocol to calculate 
the model parameters based on the experimentally de- 
termined front velocity, v(Cq)- Wc emphasize that the 
most notable feature of the model is its ability to distin- 
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guish between permanent traps (straining) and the traps 
with small but finite detachment rate. It is the latter 
traps that determine the long-time asymptotics of the 
washout curves. 

The suggested model is applicable to a wide range 
of problems in which macromolccules, stable emulsion 
drops, or pathogenic micro-organisms such as bacteria 
and viruses are transported in flow through a porous 
medium. While the model is purely phenomenological 
in nature, the mapping of the parameters with the ex- 
perimental data as a function of flow velocity and colloid 
size will shed light on the nature of trapping for particu- 
lar colloids. The model can also be extended to account 
for variations in attachment and detachment rates for 
various colloids as needed to explain the steep deposition 
profiles near the inlet of filters pj. 
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APPENDIX: VELOCITY OF AN 
INTERMEDIATE FRONT. 

Here we derive an inequality for the velocity vab of 
an intermediate front interpolating between free-particle 
concentrations Ca and Cb [Fig. p|. 



We first write the expressions for the filtering front 
velocities in clean filter, with the inlet concentrations Ca 
and C B > C A [cf. Eq. ©], 



The velocity vab of the filtering front interpolating be- 
tween Ca and Cb [Fig. p| is given by 



— -1)Ca = J2Nini(C A ), 

— -i)c B = ViWCs). 



( v 



-l)(C B - Ca) = J2 N i[ n i( C B) - ni(C A )]. 

(A.l) 



\VAB 

Combining these equations, we obtain 
C B Ca Cb — Ca 



vb v a 



VAB 



(A.2) 



Fro m he re we conclude that the left-hand side (lhs) of 
Eq. (A.2) is positive. Solving for vab and expressing the 
difference vab — V A , we have 



vab - va 



C b (vb - va) 
Cb Ca\ 



(A.3) 



. v B v A L-^ 
For the model with saturation [Eq. (p2|)1, we saw that 

vb > va, thus vab > va- 
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